home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Libris Britannia 4
/
science library(b).zip
/
science library(b)
/
MATHEMAT
/
1769.ZIP
/
PFSAF.FOR
< prev
next >
Wrap
Text File
|
1988-10-05
|
54KB
|
2,750 lines
$include: 'pfcom.for'
open(unit=7,file='algin',status='old')
open(unit=8,file='algout',status='unknown')
call datim
tbeg=tnow
do 4 i=1,nvar
4 ittb(i)=0
ios(1)=1
iosi(1)=1
ipv(1)=1
pv(1,1)='?'
ittb(1)=1
itte(1)=1
itb(1)=1
ite(1)=1
ivn(1)=2
ivp(1)=1
nb=7
nd=0
nv=1
ns=1
np=2
nt=1
maxe=9999
maxnt=0
maxite=0
itetop=1
ndump=0
minns=9999
lintot=0
idflg=1
call defh
ma=ntot
mb=nvar
mc=nterm
md=ninch
me=ninl
mg=nder
mh=nlab
mi=ndind
mj=mitay
write(8,1)ma,mb,mc,md,me,mg,mh,mi,mj
1 format(' version 85/6/21 ',
+ ' NTOT=',i6,' NVAR=',i4,' NTERM=',i5,' NINCH=',i5,/
+ ,' NINL=',i6,' NDER=',i3,' NLAB=',i3,' NDIND=',i2
+ ,' MITAY=',i3)
call monitr
call datim
write(8,90)nd,nv,ns,maxnt,maxite,lintot
90 format(' nu. derivatives=',i3,' nu. variables=',i3,
+ ' nu. sums=',i3,/,' max nu. terms=',i5,
+ ' max total nu. variables=',i6,' input lines=',i5)
t=tnow-tbeg
write(8,91)t
91 format(' time used=',f9.2,' seconds')
stop
end
subroutine add(is,js,t,ians)
$include: 'pfcom.for'
nt1=nt+1
ib=ittb(is)
ie=itte(is)
jb=ittb(js)
je=itte(js)
i=itb(ie)
if(ivn(i).ne.2 .or. ivp(i).le.maxe)goto 3
do 1 it=ib,ie
i=itb(it)
if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 2
1 continue
2 ie=it-1
3 i=itb(je)
if(ivn(i).ne.2 .or.ivp(i).le.maxe)goto 7
do 4 it=jb,je
i=itb(it)
if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 5
4 continue
5 je=it-1
7 if(ib.gt.ie)goto 70
if(jb.gt.je)goto 12
if(icomp(ib,jb))8,60,30
8 jl=jb
if(jb.eq.je)goto 10
if(icomp(ib,je))10,16,18
10 call cpynt(jb,je,t)
12 call cpynt(ib,ie,1.d0)
14 call defsum(ians,nt1,nt)
return
16 save=tm(je)
tm(je)=tm(ib)/t+tm(je)
call cpynt(jb,je,t)
tm(je)=save
ib=ib+1
jb=je+1
goto 7
18 jr=je
20 if(jr-jl.eq.1)goto 28
jc=(jl+jr)/2
if(icomp(ib,jc))22,26,24
22 jl=jc
goto 20
24 jr=jc
goto 20
26 save=tm(jc)
tm(jc)=tm(ib)/t+tm(jc)
call cpynt(jb,jc,t)
tm(jc)=save
ib=ib+1
jb=jc+1
goto 7
28 call cpynt(jb,jl,t)
jb=jr
30 il=ib
if(ib.eq.ie)goto 31
if(icomp(jb,ie))31,36,38
31 call cpynt(ib,ie,1.d0)
32 call cpynt(jb,je,t)
goto 14
36 save=tm(ie)
tm(ie)=tm(ie)+t*tm(jb)
call cpynt(ib,ie,1.d0)
tm(ie)=save
jb=jb+1
ib=ie+1
goto 70
38 ir=ie
40 if(ir-il.eq.1)goto 48
ic=(il+ir)/2
if(icomp(jb,ic))42,46,44
42 il=ic
goto 40
44 ir=ic
goto 40
46 save=tm(ic)
tm(ic)=tm(ic)+t*tm(jb)
call cpynt(ib,ic,1.d0)
tm(ic)=save
jb=jb+1
ib=ic+1
goto 7
48 call cpynt(ib,il,1.d0)
ib=ir
goto 8
60 save=tm(ib)
tm(ib)=tm(ib)+t*tm(jb)
call cpynt(ib,ib,1.d0)
tm(ib)=save
ib=ib+1
jb=jb+1
goto 7
70 if(jb.gt.je)goto 14
goto 32
end
subroutine bomb(i,h)
$include: 'pfcom.for'
character*8 h
goto(100,100,30,40,50),i
30 write(8,31)h
31 format(a8)
goto 100
40 write(8,41)h,h
41 format(a8,' is too small for your problem.',
+ ' Recompile with bigger ',a8)
goto 100
50 continue
100 write(8,101)line,ap
101 format('input line number',i5,' was:',/,127a1)
call datim
write(8,90)nd,nv,ns,maxnt,maxite,line
90 format(/,'nu. derivatives=',i3,
+ ' nu. variables=',i3,'nu. sums=',i3,' max nu. terms=',
+ i5,/,' max total nu. variables=',i6,' input lines=',i5)
if(ndump.ne.0)call dump
t=tnow-tbeg
write(8,91)t
91 format('time used=',f9.2,' seconds')
stop
end
subroutine bkup(j,k)
$include: 'pfcom.for'
if(j.eq.k)goto 2
do 1 i=j,k-1
itb(i)=itb(i+1)
ite(i)=ite(i+1)
1 tm(i)=tm(i+1)
2 k=k-1
return
end
subroutine call
$include: 'pfcom.for'
call nxt(3,nchpl,' ',j1)
call numl(j1,nchpl,j2,j3,j4)
if(linlab(j4).eq.-1)call bomb(3,'call ')
linret=line
line=linlab(j4)
return
end
subroutine coef(is,j3,c)
c c=coefficient, j3=start of variable (on return)
$include: 'pfcom.for'
f=1.d0
sign=1.d0
call nxtnb(is,nchpl,j2)
call srch(j2,nchpl,'*',jj)
if(jj.eq.-1)jj=nchpl
call srch(j2,nchpl,'/',jk)
if(jk.eq.-1)jk=nchpl
call srchch(j2,nchpl,j3)
j3=min0(jj,jk,j3)
if(a(j2).eq.'-')goto 4
if(a(j2).eq.'+')goto 10
j4=j2
goto 12
4 sign=-1.d0
10 call nxtnb(j2+1,nchpl,j4)
12 c=sign
if(a(j4).eq.'.')goto 15
call intval(j4,i)
if(i.eq.-1 .and. j4.ne.j3)call bomb(3,'coef1 ')
if(i.eq.-1)return
15 c=0.d0
do 30 j1=j4,nchpl
if(a(j1).eq.'.' .and. f.lt..5d0)call bomb(3,'coef2 ')
if(a(j1).eq.'.')f=.1d0
if(a(j1).eq.'.')goto 30
call intval(j1,i)
if(i.eq.-1)goto 40
if(f.gt..5d0)c=10.d0*c+i
if(f.lt..5d0)c=c+i*f
if(f.lt..5d0)f=f*.1d0
30 continue
40 c=c*sign
return
end
subroutine comput
$include: 'pfcom.for'
300 call nxt(1,nchpl,' ',j1)
call numv(j1,nchpl,j2,j3,j4)
call nxt(j3,nchpl,'=',j5)
call srch(j5,nchpl,'/',j6)
if(j6.eq.-1)goto 320
c /
call nxtnb(j5+1,j6-1,j7)
if(a(1).eq.'d')goto 301
call numv(j7,j6-1,j8,j9,j10)
call numv(j6+1,nchpl,j11,j12,j13)
if(ittb(j10).eq.0 .or. ittb(j13).eq.0)call bomb(3,'not sum ')
call recip(j13,1)
nt1=nt+1
call mult(ittb(j10),itte(j10),ittb(1),itte(1))
call defsum(j4,nt1,nt)
call outqu(j4)
minns=min0(minns,ios(1))
it=ittb(1)
itte(1)=it
ite(it)=itb(it)
k=itb(it)
ivn(k)=2
ivp(k)=0
call pack
return
c compute sj4=dsj10/dfj14
301 call numv(j7+1,j6-1,j8,j9,j10)
if(a(j7).ne.'d')call bomb(3,'comput1 ')
if(idflg.ne.0)call dertab
call nxt(j6,nchpl,'d',j11)
call numv(j11+1,nchpl,j12,j13,j14)
call srch(j13,nchpl,'=',j15)
iord=1
if(j15.eq.-1)goto 311
call nxtnb(j15+1,nchpl,j16)
call intval(j16,j17)
if(j17.ne.-1)goto 360
call numv(j16,nchpl,j19,j20,j18)
if(ittb(j18).eq.0)call bomb(3,'comput14')
iord=1.000000001d0*tm(ittb(j18))
if(iord.gt.0)goto 311
call bomb(3,'comput15')
360 call getexp(j15+1,iord,j16,0)
if(j16.ne.-1 .or. iord.le.0)call bomb(3,'comput2 ')
311 do 313 ii=1,iord
if(ii.eq.1)call deriv(j10,j14,1)
if(ii.ne.1)call deriv(1,j14,1)
call outqu(1)
if(tm(ittb(1)).eq.0.)goto 314
313 continue
314 call swap(j4)
return
c not /
320 call srch(j5,nchpl,'+',j6)
if(a(1).eq.'d')call bomb(3,'comput3 ')
if(j6.ne.-1)goto 321
call srch(j5,nchpl,'-',j6)
tmm=-1.d0
if(j6.ne.-1)goto 322
call srch(j5,nchpl,'*',j6)
if(j6.ne.-1)goto 339
call srch(j5,nchpl,'^',j6)
if(j6.ne.-1)goto 339
c compute sj4=sj12
call numv(j5,nchpl,j6,j11,j12)
if(ittb(j12).eq.0)call bomb(3,'comput4 ')
call copys(j12,j4)
call outqu(j4)
return
c compute sj4=sj9 + or - sj12
321 tmm=1.d0
322 call numv(j5+1,j6-1,j7,j8,j9)
call numv(j6+1,nchpl,j10,j11,j12)
call add(j9,j12,tmm,j4)
call outqu(j4)
return
c compute sj4=sj9*sj12
339 call numv(j5+1,j6-1,j7,j8,j9)
if(a(j6).eq.'^')goto 350
if(a(j6+1).eq.'*')goto 349
call numv(j6+1,nchpl,j10,j11,j12)
nt1=nt+1
call mult(ittb(j9),itte(j9),ittb(j12),itte(j12))
call defsum(j4,nt1,nt)
call outqu(j4)
return
349 j6=j6+1
350 call numv(j6+1,nchpl,j10,j11,j12)
call expnd(j9,j12,j4)
return
end
subroutine copy(a,b,n)
character*1 a,b
dimension a(2),b(2)
do 1 i=1,n
1 b(i)=a(i)
return
end
subroutine copys(j1,j2)
$include: 'pfcom.for'
nt1=nt+1
call cpynt(ittb(j1),itte(j1),1.d0)
call defsum(j2,nt1,nt)
return
end
subroutine count
$include: 'pfcom.for'
call nxt(3,nchpl,' ',j1)
call numv(j1,nchpl,j3,j4,j5)
do 12 i=nchpl,j4+1,-1
if(a(i).ne.' ')goto 14
12 continue
call bomb(3,'count1 ')
14 do 16 j=i,j4+1,-1
if(a(j).eq.' ')goto 18
16 continue
18 call numv(j,i,j6,j7,j8)
if(ittb(j8).eq.0)call bomb(3,'count2 ')
if(ittb(j5).ne.0)goto 20
k=itetop+1
call incnt
itb(nt)=k
call nite(k)
tm(nt)=itte(j8)+1-ittb(j8)
ivn(k)=2
ivp(k)=0
call defsum(j5,nt,nt)
return
20 ia=ittb(j5)
ib=itte(j5)
k=itb(ia)
l=ite(ia)
ivn(k)=2
ivp(k)=0
itte(j5)=ia
ite(ia)=k
tm(ia)=itte(j8)+1-ittb(j8)
if(ia.ne.ib .or. k.ne.l)minns=min0(minns,ios(j5))
call pack
return
end
subroutine cpynt(ia,ib,tmm)
c copy terms ia-ib to nt+1,..., coefficients are multiplied by tmm
$include: 'pfcom.for'
if(ib.lt.ia)return
if(ia.lt.1 .or. ib.gt.nt)call bomb(3,'cpynt2 ')
do 15 i=ia,ib
ja=itb(i)
if(ivn(ja).eq.2 .and. ivp(ja).gt.maxe)goto 15
k=itetop
call incnt
tm(nt)=tm(i)*tmm
itb(nt)=k+1
do 10 j=ja,ite(i)
k=k+1
ivn(k)=ivn(j)
10 ivp(k)=ivp(j)
call nite(k)
15 continue
16 return
end
subroutine datim
$include: 'pfcom.for'
call date(iy,mon,id)
call time(ih,min,is,if)
tnow=is+60.*(min+60.*ih)+if/100.
write(8,1)iy,mon,id,ih,min,is
1 format(' job run ',i4,'/',i2,'/',i2,i9,':',i2,':',i2)
return
end
subroutine datime
$include: 'pfcom.for'
call time(ih,min,is,if)
tnow=is+60.*(min+60.*ih)+if/100.
return
end
subroutine defc(a,n)
character*1 a(2)
do 1 i=2,n
j=a(i-1)
1 a(i)=j+1
return
end
subroutine defh
$include: 'pfcom.for'
character*1 ga(26),gc(26),gn(10)
data ga/'a','b','c','d','e','f','g','h','i','j','k','l','m'
, ,'n','o','p','q','r','s','t','u','v','w','x','y','z'/
data gc/'A','B','C','D','E','F','G','H','I','J','K','L','M'
, ,'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
data gn/'0','1','2','3','4','5','6','7','8','9'/
do 10 i=1,10
10 hn(i)=gn(i)
do 20 i=1,26
ha(i)=ga(i)
20 hc(i)=gc(i)
return
end
subroutine defsum(is,ittbv,ittev)
$include: 'pfcom.for'
if(ittev.ne.nt .or. ittev+1.ne.ittbv)goto 1
call niltrm
ittev=nt
1 if(ittbv.gt.ittev)call bomb(3,'defsum1 ')
if(ittb(is).eq.0)goto 20
j=ios(is)
minns=min0(minns,j)
if(ittbv.ge.ittb(is) .and. ittev.le.itte(is))goto 8
if(j.eq.ns .and. ittbv.ge.ittb(is))goto 8
if(ittbv.gt.itte(iosi(ns)))goto 4
call bomb(3,'defsum2 ')
4 ios(is)=ns
do 6 i=j,ns-1
k=iosi(i+1)
iosi(i)=k
6 ios(k)=i
iosi(ns)=is
8 ittb(is)=ittbv
itte(is)=ittev
call pack
return
20 if(ittbv.le.itte(iosi(ns)))call bomb(3,'defsum3 ')
ns=ns+1
ios(is)=ns
iosi(ns)=is
ittb(is)=ittbv
itte(is)=ittev
minns=min0(minns,ns)
call pack
return
end
subroutine delet
$include: 'pfcom.for'
c delete sj4
call nxt(4,nchpl,' ',j1)
1 call numv(j1,nchpl,j2,j3,j4)
call delsum(j4)
if(j3.eq.nchpl)goto 10
call srchnb(j3+1,nchpl,j1)
if(j1.ne.-1)goto 1
10 return
end
subroutine delsum(is)
$include: 'pfcom.for'
if(ittb(is).eq.0)return
j=ios(is)
minns=min0(minns,j)
if(j.eq.0 .or. ns.eq.0)call bomb(3,'delsum1 ')
ns=ns-1
if(j.gt.ns)goto 10
do 2 i=j,ns
k=iosi(i+1)
iosi(i)=k
2 ios(k)=i
10 iosi(ns+1)=0
ios(is)=0
ittb(is)=0
itte(is)=0
call pack
return
end
subroutine deriv(j10,ind,jans)
c d(sum j10)/d(var ind) is put in sum j4
c ider(1,nd) is dependent variable in derivative
c ider(2,nd) is independent variable in derivative
c ider(3,nd) is variable equal to derivative
$include: 'pfcom.for'
character*1 bl
common/i/if
nt1=nt+1
if(ittb(j10).lt.1 .or. itte(j10).gt.nt)call bomb(3,'deriv2 ')
do 90 it=ittb(j10),itte(j10)
ia=itb(it)
ib=ite(it)
do 85 ic=ia,ib
ie=ivn(ic)
if(ie.ne.ind)goto 40
call incnt
itb(nt)=itetop+1
tm(nt)=tm(it)*ivp(ic)
if=itetop
do 36 i=ia,ib
if(i.eq.ic)goto 32
if=if+1
ivn(if)=ivn(i)
ivp(if)=ivp(i)
goto 36
32 if(ivp(i).eq.1 .and. ia.ne.ib)goto 36
if=if+1
ivn(if)=ivn(i)
ivp(if)=ivp(i)-1
36 continue
call nite(if)
if(if.lt.itb(nt))call bomb(3,'deriv3 ')
goto 85
40 if(nd.eq.0)goto 85
i1=1
41 j1=idt(i1,ie)
if(j1.eq.0)goto 85
k1=ider(2,j1)
level=1
if(k1.eq.ind)goto 50
i2=1
42 j2=idt(i2,k1)
if(j2.eq.0)goto 77
k2=ider(2,j2)
level=2
if(k2.eq.ind)goto 50
i3=1
43 j3=idt(i3,k2)
if(j3.eq.0)goto 76
k3=ider(2,j3)
level=3
if(k3.eq.ind)goto 50
i4=1
44 j4=idt(i4,k3)
if(j4.eq.0)goto 75
k4=ider(2,j4)
level=4
if(k4.eq.ind)goto 50
i5=1
45 j5=idt(i5,k4)
if(j5.eq.0)goto 74
k5=ider(2,j5)
level=5
if(k5.eq.ind)goto 50
l1=ipv(ie)
l2=ipv(k1)
l3=ipv(k2)
l4=ipv(k3)
l5=ipv(k4)
l6=ipv(k5)
bl=' '
write(8,49)(pv(i,ie),i=1,l1),bl,(pv(i,k1),i=1,l2),bl
+ ,(pv(i,k2),i=1,l3),bl,(pv(i,k3),i=1,l4),bl,(pv(i,k4),i=1,l5)
+ ,bl,(pv(i,k5),i=1,l6)
49 format('chain too long: ',90a1)
call bomb(3,'deriv4 ')
50 call incnt
if=itetop
itb(nt)=if+1
tm(nt)=tm(it)*ivp(ic)
do 55 i=ia,ib
if(i.eq.ic)goto 52
if=if+1
ivn(if)=ivn(i)
ivp(if)=ivp(i)
goto 55
52 if(ivp(i).eq.1)goto 55
if=if+1
ivn(if)=ivn(i)
ivp(if)=ivp(i)-1
55 continue
call deriv2(j1)
if(level.eq.1)goto 77
call deriv2(j2)
if(level.eq.2)goto 76
call deriv2(j3)
if(level.eq.3)goto 75
call deriv2(j4)
if(level.eq.4)goto 74
call deriv2(j5)
i5=i5+1
if(i5.le.ndind)goto 45
74 i4=i4+1
if(i4.le.ndind)goto 44
75 i3=i3+1
if(i3.le.ndind)goto 43
76 i2=i2+1
if(i2.le.ndind)goto 42
77 i1=i1+1
if(i1.le.ndind)goto 41
85 continue
90 continue
call defsum(jans,nt1,nt)
call order(jans)
return
end
subroutine deriv2(j)
$include: 'pfcom.for'
common/i/if
if=if+1
call nite(if)
ivn(if)=ider(3,j)
ivp(if)=1
return
end
subroutine dertab
$include: 'pfcom.for'
idflg=0
do 302 i=1,ndind
do 302 j=1,nvar
302 idt(i,j)=0
do 304 i=1,nd
idepv=ider(1,i)
do 303 j=1,ndind
if(idt(j,idepv).eq.0)goto 304
303 continue
call bomb(4,'ndind ')
304 idt(j,idepv)=i
return
end
subroutine dfine
$include: 'pfcom.for'
call nxt(4,nchpl,' ',j1)
call nxtch(j1,nchpl,j2)
call numv(j2,nchpl,j4,j5,j6)
ittb6=nt+1
call srch(j5,nchpl,'=',j7)
if(j7.eq.-1)goto 410
call term(j7+1)
goto 412
410 call read(ib)
if(np.ge.2)call prntap
if(ib.eq.1)goto 412
call term(1)
goto 410
412 call defsum(1,ittb6,nt)
call order(1)
call outqu(1)
call swap(j6)
return
end
subroutine dump
$include: 'pfcom.for'
character*1 bl
data bl/' '/
write(8,1)nd,nv,ns,np,nt,line,maxe,maxnt,maxite
1 format(/,'nd=',i3,' nv=',i3,' ns=',i3,' np=',i2,' nt=',i5,
+ ' line=',i3,' maxe=',i2,' maxnt=',i5,' maxite=',i6)
do 10 i=1,nv
j=ipv(i)
if(j.eq.nch)goto 3
jp=j+1
do 4 jj=jp,nch
4 pv(jj,i)=bl
3 write(8,5)i,(pv(k,i),k=1,nch),ittb(i),itte(i),ios(i),iosi(i)
5 format('i=',i5,' pv=',15a1,' ittb=',i5,' itte=',i5
+ ,' ios=',i3,' iosi=',i3)
10 continue
write(8,11)
11 format(' ')
write(8,13)
13 format(' i tm itb ite')
do 20 i=1,nt,4
j=i+1
k=i+2
l=i+3
write(8,21)i,tm(i),itb(i),ite(i),j,tm(j),itb(j),ite(j)
+ ,k,tm(k),itb(k),ite(k),l,tm(l),itb(l),ite(l)
21 format(4(i9,f8.2,2i6))
20 continue
write(8,11)
do 22 i=1,nd
j=ider(3,i)
l=ipv(j)
22 write(8,24)ider(1,i),ider(2,i),(pv(m,j),m=1,l)
24 format('the derivative of',i4,' wrt',i4,' is ',15a1)
write(8,11)
write(8,40)(i,ivn(i),ivp(i),i=1,itetop)
40 format(10(1x,i5,2i3))
return
end
subroutine expnd(jb,jp,jans)
$include: 'pfcom.for'
ia=ittb(jb)
if(ia.eq.0)call bomb(3,'expnd1 ')
if(tm(ia).ne.1.d0)call bomb(3,'expnd2 ')
ic=itb(ia)
if(ivn(ic).ne.2 .or. ivp(ic).ne.0)call bomb(3,'expnd3 ')
ic=ia+1
id=itte(jb)
5 call cpynt(ia,ia,1.d0)
if(ic.le.id)goto 10
6 call defsum(jans,nt,nt)
return
10 ig=nt
ie=ittb(jp)
tp=tm(ie)
if(tp.eq.0.d0)goto 6
jprvb=ig
jprve=ig
tprv=1.d0
do 20 kp=1,nb
if(dabs(tp).lt.1.d-11)goto 25
tprv=tprv*tp/kp
tmt(kp)=tprv
tp=tp-1.d0
itbt(kp)=nt+1
call mult(jprvb,jprve,ic,id)
itet(kp)=nt
if(itbt(kp).gt.nt)goto 25
jprvb=itbt(kp)
20 jprve=nt
kp=nb+1
25 kp=kp-1
do 40 i=1,kp
ih=itbt(i)
ii=itet(i)
tp=tmt(i)
do 30 il=ih,ii
30 tm(il)=tm(il)*tp
40 continue
call defsum(jans,ig,ii)
call order(jans)
return
end
subroutine extrct
$include: 'pfcom.for'
call nxt(2,nchpl,' ',j1)
call srch(j1,nchpl,'-',j6)
if(j6.ne.-1)goto 300
num=1
do 1 i=2,10
if(a(j1-1).eq.hn(i))num=i-1
1 continue
call numv(j1,nchpl,j3,j4,j5)
call nxt(j4,nchpl,'=',j6)
call nxt(j6,nchpl,'m',j7)
call numv(j7+1,nchpl,j8,j9,j10)
k=ittb(j10)
if(k.eq.0)call bomb(3,'extrct1 ')
kk=tm(k)+.1
if(kk.gt.0)goto 15
12 write(8,14)tm(k),(a(i),i=j8,j9)
14 format('bad term number ',d12.3,5x,15a1)
call bomb(3,'extrct2 ')
15 call nxt(j9+1,nchpl,'o',j11)
call numv(j11+2,nchpl,j12,j13,j14)
kl=itte(j14)+1-ittb(j14)
if(kk.gt.kl)goto 12
km=ittb(j14)+kk-1
nt1=nt+1
ktop=min0(km+num-1,itte(j14))
call cpynt(km,ktop,1.d0)
call defsum(j5,nt1,nt)
call outqu(j5)
return
300 call numv(j1,nchpl,j2,j3,j4)
call nxt(j3,j6,'=',j5)
call nxt(j5,j6,'t',j9)
call nxt(j9,j6,' ',j10)
call nxtnb(j10,j6,j11)
call srch(j11,j6,' ',j12)
if(j12.eq.-1)j12=j6
call getint(j11,j12-1,j13)
call srchnb(j6+1,nchpl,j8)
call srch(j8,nchpl,' ',j14)
call getint(j8,j14-1,j15)
call nxt(j14,nchpl,'o',j19)
if(a(j19+1).ne.'f' .and. a(j19+2).ne.' ')
+ call bomb(3,'extrct5 ')
330 call numv(j19+2,nchpl,j16,j17,j18)
if(ittb(j18).eq.0)call bomb(3,'extrct6')
if(j15.gt.itte(j18)-ittb(j18)+1)j15=itte(j18)-ittb(j18)+1
if(j15.lt.j13)call bomb(3,'extrct7 ')
nt1=nt+1
call cpynt(ittb(j18)+j13-1,ittb(j18)+j15-1,1.d0)
call defsum(j4,nt1,nt)
call outqu(j4)
return
end
subroutine fmtc(it,if,il,posch)
$include: 'pfcom.for'
character*1 aa,pa,posch
common/pr1/pa(232),aa(22)
character*127 cdum
common/fmt/cdum
il=22
tt=tm(it)
if(tt.ne.0.)goto 3
aa(9)='0'
if=9
il=9
return
3 ab=tm(it)*64.d0
if(ab.gt.0.d0)j=ab+.5d0
if(ab.lt.0.d0)j=ab-.5d0
t=j/64.d0
ac=dmax1(1.d0,dabs(tm(it)))
if(dabs((t-tt)/ac).lt.1.d-11)tm(it)=t
tt=tm(it)
if(dabs(tt).gt.9999.9d0)goto 16
6 format(22a1)
14 if(dabs(tt).gt.1.d0)goto 20
if(dabs(tt).gt..001d0)goto 24
16 write(cdum,18)tt
18 format(d22.13)
read(cdum,6)aa
ig=22
181 ig=ig-1
if(aa(ig).ne.'E' .and. aa(ig).ne.'e')goto 181
ih=ig
182 ih=ih-1
if(aa(ih).eq.'0')goto 182
do 183 i=ig,22
183 aa(ih+1+i-ig)=aa(i)
il=ih+1+22-ig
goto 32
20 write(cdum,22)tt
22 format(f22.12)
goto 28
24 write(cdum,26)tt
26 format(f22.16)
28 read(cdum,6)aa
30 if(aa(il).ne.'0')goto 32
il=il-1
goto 30
32 do 34 if=1,il
if(aa(if).ne.' ')goto 36
34 continue
36 if(tt.gt.0.d0)if=if-1
if(tt.gt.0.d0)aa(if)=posch
return
end
subroutine fort(it,ib)
$include: 'pfcom.for'
character*1 aa,pa
common/pr1/pa(232),aa(22)
character*127 cdum
common/fmt/cdum
if(np.eq.0)return
j=ib
call fmtc(it,if,il,'+')
do 42 i=if,il
ib=ib+1
42 a(ib)=aa(i)
ic=itb(it)
id=ite(it)
do 60 iv=ic,id
ip=iabs(ivp(iv))
in=ivn(iv)
ie=ipv(in)
iq=0
if(ip.eq.1)goto 47
write(cdum,43)ip
43 format(i20)
read(cdum,6)aa
6 format(22a1)
il=20
do 44 if=1,20
if(aa(if).ne.' ')goto 45
44 continue
45 iq=il+1-if+2
47 if(iq+ie+1.le.72)goto 52
j=0
write(8,50)(a(ich),ich=1,ib)
50 format(80a1)
do 51 i=1,80
51 a(i)=' '
a(6)='.'
ib=11
52 ib=ib+1
a(ib)='*'
if(ivp(iv).lt.0)a(ib)='/'
do 53 i=1,ie
ib=ib+1
53 a(ib)=pv(i,in)
if(iq.eq.0)goto 60
ib=ib+1
a(ib)='*'
ib=ib+1
a(ib)='*'
do 54 i=if,il
ib=ib+1
54 a(ib)=aa(i)
60 continue
k=0
if(a(j+2).eq.'1' .and. a(j+3).eq.'.' .and. a(j+4).eq.'*')k=1
if(k.eq.0)goto 61
a(j+2)=' '
a(j+3)=a(j+1)
a(j+4)=' '
a(j+1)=' '
61 write(8,50)(a(ich),ich=1,ib)
return
end
subroutine fortrn
$include: 'pfcom.for'
call nxt(1,nchpl,' ',j1)
call numv(j1+1,nchpl,j2,j3,j4)
if(np.eq.0 .or. ittb(j4).eq.0)goto 10
do 1 i=1,nchpl
1 a(i)=' '
ia=ipv(j4)
ib=6
do 2 i=1,ia
ib=ib+1
2 a(ib)=pv(i,j4)
ib=ib+1
a(ib)='='
ic=ittb(j4)
id=itte(j4)
j=0
do 4 it=ic,id
call fort(it,ib)
do 3 i=1,nchpl
3 a(i)=' '
a(6)='.'
ib=7
j=j+1
if(j.lt.10)goto 4
a(6)=' '
call copy(pv(1,j4),a(7),ia)
a(7+ia)='='
call copy(pv(1,j4),a(8+ia),ia)
ib=7+ia+ia
j=0
4 continue
10 return
end
subroutine getint(j1,j2,j3)
c on return j3=integer value of a(j1)-a(j2)
j3=0
if(j2.lt.j1)call bomb(3,'getint1 ')
do 1 i=j1,j2
call intval(i,j)
1 j3=10*j3+j
return
end
subroutine getexp(j1,ie,k,if)
c begin looking at a(j1), exponent put in ie, a(k) is next nonblank
c (or k=-1 if only blanks are found)
$include: 'pfcom.for'
ie=1
call srchnb(j1,nchpl,k)
if(k.eq.-1)return
if(a(k).eq.'*'.and. a(k+1).ne.'*')return
ise=1
if(a(k).eq.'*'.and. a(k+1).eq.'*')goto 35
if(a(k).eq.'/')return
call srchch(j1,nchpl,j2)
if(j2.eq.k .and. if.ne.0)return
if(j2.ne.k)goto 9
32 call numv(j2,nchpl,j3,j4,j5)
if(ittb(j5).eq.0)call bomb(3,'getexp10')
tmittb=tm(ittb(j5))
ie=tmittb+.5d0
if(tmittb.lt.0.d0)ie=tmittb-.5d0
re=ie
if(dabs(re-tmittb).gt.1.d-11)call bomb(3,'getexp11')
k=-1
if(j4.ge.nchpl-1)return
call srchnb(j4+1,nchpl,k)
return
9 if(a(k).eq.'+')goto 20
if(a(k).eq.'-')goto 10
if(a(k).eq.'^')goto 30
j3=k
1 me=0
5 call intval(j3,j4)
me=10*me+j4
j3=j3+1
if(j3.gt.nchpl)goto 70
if(a(j3).eq.' ')goto 50
if(a(j3).eq.'*')goto 60
if(a(j3).eq.'/')goto 60
goto 5
10 ise=-1
20 if(k+1.gt.nchpl)call bomb(3,'getexp2 ')
call nxtnb(k+1,nchpl,j3)
if(j3.eq.-1)call bomb(3,'getexp4 ')
goto 1
30 call nxtnb(k+1,nchpl,j3)
31 k=j3
j2=j3
call srchch(j3,nchpl,j4)
if(j3.eq.j4)goto 32
if(a(k).eq.'-')goto 10
if(a(k).eq.'+')goto 20
goto 1
35 if(k.ge.nchpl-1)call bomb(3,'getexp6 ')
call nxtnb(k+2,nchpl,j3)
if(j3.eq.-1)call bomb(3,'getexp8 ')
goto 31
50 ie=me*ise
call srchnb(j3,nchpl,k)
return
60 k=j3
if(j3.eq.nchpl)call bomb(3,'getexp9 ')
ie=me*ise
return
70 ie=me*ise
k=-1
return
end
function icomp(ia,ib)
$include: 'pfcom.for'
c if ia should precede ib icomp=1
c if ib should precede ia icomp=-1
c if terms are the same except possibly coefficients, icomp=0
ja=itb(ia)
jb=itb(ib)
ka=ite(ia)
kb=ite(ib)
1 if(ivn(ja)-ivn(jb))2,8,9
2 if(ivp(ja))3,10,4
4 icomp=-1
return
5 if(jb.gt.kb)goto 7
3 icomp=1
return
6 ja=ja+1
11 jb=jb+1
if(ja.gt.ka)goto 5
if(jb.gt.kb)goto 4
goto 1
7 icomp=0
return
8 if(ivp(ja)-ivp(jb))3,6,4
9 if(ivp(jb))4,11,3
10 ja=ja+1
if(ja.gt.ka)goto 3
goto 1
end
subroutine iftest(go)
$include: 'pfcom.for'
character*1 aj1,aj2
go=0.
call nxt(2,nchpl,' ',j1)
call nxt(j1,nchpl,'.',j)
call numv(j1,j-1,j3,j4,j5)
aj1=a(j+1)
aj2=a(j+2)
if(aj1.eq.'n')goto 20
if(aj1.eq.'l' .or. aj1.eq.'g')goto 22
if(aj1.ne.'e')call bomb(3,'iftest1 ')
if(aj2.ne.'q')call bomb(3,'iftest2 ')
4 if(a(j+3).ne.'.')call bomb(3,'iftest3 ')
call numv(j+4,nchpl,j6,j7,j8)
if(j5.ne.j8)goto 10
if(aj1.eq.'n')return
5 call nxtch(j7+1,nchpl,j9)
do 6 j=j9,nchpl
a(j+1-j9)=a(j)
6 a(j)=' '
go=1.d0
call copy(a,ap,nchpl)
return
10 k1=itb(ittb(j5))
k2=ite(itte(j5))
k3=itb(ittb(j8))
k4=ite(itte(j8))
t5=tm(ittb(j5))
t8=tm(ittb(j8))
if(aj1.eq.'l')goto 23
if(aj1.eq.'g')goto 25
if(k2-k1.ne.k4-k3)goto 30
k=k3
do 11 i=k1,k2
if(ivn(i).ne.ivn(k))goto 30
if(ivp(i).ne.ivp(k))goto 30
11 k=k+1
k1=ittb(j5)
k2=itte(j5)
k3=ittb(j8)
k4=itte(j8)
if(k2-k1.ne.k4-k3)goto 30
k=k3
do 12 i=k1,k2
den=dabs(tm(i))+dabs(tm(k))
if(den.lt.1.d-11)goto 12
if(dabs(tm(i)-tm(k))/den.gt.1.d-11)goto 30
12 k=k+1
if(aj1.eq.'n')return
goto 5
20 if(a(j+2).ne.'e')call bomb(3,'iftest4 ')
goto 4
22 if(aj2.eq.'e' .or. aj2.eq.'t')goto 4
call bomb(5,'iftest5 ')
23 if(aj2.eq.'e')goto 24
if(t5.lt.t8)goto 5
return
24 if(t5.le.t8)goto 5
return
25 if(aj2.eq.'e')goto 26
if(t5.gt.t8)goto 5
return
26 if(t5.ge.t8)goto 5
return
30 if(aj1.eq.'n')goto 5
return
end
subroutine incnt
$include: 'pfcom.for'
nt=nt+1
maxnt=max0(nt,maxnt)
if(maxnt.gt.nterm)call bomb(4,'nterm ')
return
end
subroutine intval(j1,i)
c on return i=integer value of a(j1)
$include: 'pfcom.for'
do 1 j=1,10
if(a(j1).eq.hn(j))goto 2
1 continue
i=-1
return
2 i=j-1
return
end
subroutine monitr
$include: 'pfcom.for'
call readin
ntp=0
tprev=tnow
10 if(ha(1).ne.'a')call bomb(3,'monitr1 ')
if(ittb(2).ne.0)call bomb(3,'useorder')
if(ha(1).eq.'a')goto 11
if(ns.lt.2)goto 14
do 20 ms=2,ns
is=iosi(ms-1)
js=iosi(ms)
if(is.eq.0 .or. js.eq.0)call bomb(3,'mon1 ')
if(ios(is).ne.ms-1 .or. ios(js).ne.ms)call bomb(3,'mon2 ')
if(ittb(is).eq.0 .or. ittb(js).eq.0)call bomb(3,'mon3 ')
ia=itte(is)
ib=ittb(js)
if(ite(ia)+1.ne.itb(ib))call bomb(3,'mon4 ')
do 16 it=ib,itte(js)
if(itb(it).gt.ite(it))call bomb(3,'mon5 ')
if(dabs(tm(it)).gt.1.d-11)goto 15
if(itb(it).ne.ite(it))call bomb(3,'mon6 ')
if(ivn(itb(it)).ne.2)call bomb(3,'mon7 ')
if(itb(it).gt.ite(it))call bomb(3,'mon8 ')
15 do 17 i=itb(it),ite(it)
if(ivn(i).eq.1)call bomb(3,'mon9 ')
if(ittb(ivn(i)).ne.0)call bomb(3,'mon10 ')
if(i.eq.itb(it))goto 17
if(ivp(i).eq.0)call bomb(3,'mon11 ')
if(ivn(i-1).ge.ivn(i))call bomb(3,'mon12 ')
17 continue
16 continue
20 continue
14 if(ns.eq.0)call bomb(3,'mon13 ')
ic=iosi(ns)
if(itte(ic).ne.nt)call bomb(3,'mon14 ')
if(ite(itte(ic)).ne.itetop)call bomb(3,'mon15 ')
11 call read(ib)
if(ib.eq.1 .and. np.ge.1)call prntap
if(ib.eq.1)goto 10
call datime
t=tnow-tprev
ttot=tnow-tbeg
tprev=tnow
ntmntp=nt-ntp
if(np.ge.3)write(8,3)t,ttot,ntmntp,nt,itetop,maxnt,maxite
3 format(10x,'t=',f6.2,' total t=',f8.1,' ntmntp=',i5,' nt=',i5,
+ ' itetop=',i6,' maxnt=',i5,' maxite=',i6)
if(np.ge.2)call prntap
5 ntp=nt
call copy(a,ap,nchpl)
if(a(1).eq.'c')goto 300
if(a(1).eq.'d')goto 400
if(a(1).eq.'e')goto 500
if(a(1).eq.'f')goto 600
if(a(1).eq.'g')goto 700
if(a(1).eq.'i')goto 900
if(a(1).eq.'l')goto 10
if(a(1).eq.'m')goto 1300
if(a(1).eq.'n')goto 1400
if(a(1).eq.'o')goto 1500
if(a(1).eq.'p')goto 1600
if(a(1).eq.'r')goto 1800
if(a(1).eq.'s')goto 1900
if(a(1).eq.'t')goto 2000
if(a(1).eq.'z')goto 2600
call bomb(3,'bad cmnd')
300 if(a(2).eq.'o')goto 310
if(a(2).ne.'a')call bomb(3,'monitr2 ')
call call
goto 10
310 if(a(3).eq.'m')goto 320
if(a(3).ne.'u')call bomb(3,'monitr3 ')
call count
goto 10
320 call comput
call pack
goto 10
400 if(a(2).eq.'u')goto 470
if(a(2).eq.'i')goto 320
if(a(2).ne.'e')call bomb(3,'monitr4 ')
if(a(3).eq.'l')goto 450
if(a(3).eq.'f')goto 440
call bomb(3,'monitr5 ')
440 call dfine
goto 10
450 call delet
call pack
goto 10
470 call dump
goto 10
500 if(a(2).eq.'n')return
if(a(2).ne.'x')call bomb(3,'monitr6 ')
call extrct
goto 10
600 if(a(2).eq.'l')goto 610
call fortrn
goto 10
610 call readin
goto 10
700 call nxt(1,nchpl,' ',j9)
call numl(j9,nchpl,j10,j11,j12)
if(linlab(j12).eq.-1)write(8,710)
if(linlab(j12).eq.-1)call bomb(3,'monitr7 ')
710 format('label not found')
line=linlab(j12)-1
goto 10
900 call iftest(go)
if(go.eq.1.d0)goto 5
goto 10
1300 call srch(1,nchpl,'=',j)
call srchnb(j+1,nchpl,k)
call getexp(k,maxe,j,0)
goto 10
1400 call nxt(3,nchpl,'=',j1)
call getexp(j1+1,j3,j2,0)
if(a(2).eq.'d')goto 1405
if(a(2).eq.'p')goto 1410
if(a(2).eq.'b')goto 1415
call bomb(3,'monitr8 ')
1405 if(a(3).eq.'u')goto 1406
call bomb(3,'monitr9 ')
1406 ndump=j3
goto 10
1410 np=j3
goto 10
1415 nb=j3
goto 10
1500 if(nv.ne.1)call bomb(3,'putfirst')
call ordset
goto 10
1600 if(a(3).eq.'i')goto 1620
if(a(3).ne.'o')call bomb(3,'monitr10')
1605 call read(ib)
if(np.ge.2)call prntap
if(a(1).ne.'r' .or. a(2).ne.'e' .or. a(3).ne.'t')goto 1605
goto 10
1620 call nxt(4,nchpl,' ',j1)
call numv(j1,nchpl,j2,j3,j4)
if(np.eq.1)write(8,1)
if(np.eq.1)call prntap
1 format(130a1)
call print(j4)
goto 10
1800 if(a(2).ne.'e' .or. a(3).ne.'t')call bomb(3,'return? ')
line=linret
goto 10
1900 call sub
call pack
goto 10
2000 if(a(2).eq.'a')goto 2010
call thderv
goto 10
2010 call taylor
goto 10
2600 call numv(2,nchpl,j1,j2,j3)
call zsub(j3)
call pack
goto 10
end
subroutine mult(ia,ib,ic,id)
c multiply terms ia-ib by terms ic-id
c result goes to nt+1,.....
$include: 'pfcom.for'
common/nwt/iitb,iite,jitb,jite,ntb,if
if(ia.eq.0 .or. ib.eq.0 .or. ic.eq.0 .or. id.eq.0)
+ call bomb(3,'mult1 ')
ntb=nt+1
if=0
do 20 i=ia,ib
iitb=itb(i)
iite=ite(i)
do 10 j=ic,id
jitb=itb(j)
jite=ite(j)
if(ivn(iitb)+ivn(jitb).eq.4 .and. ivp(iitb)+ivp(jitb).gt.maxe)
+ goto 20
call sub2(tm(i)*tm(j))
10 continue
20 continue
return
end
subroutine niltrm
$include: 'pfcom.for'
k=itetop+1
call incnt
tm(nt)=0.d0
itb(nt)=k
call nite(k)
ivn(k)=2
ivp(k)=0
return
end
subroutine nite(j)
$include: 'pfcom.for'
itetop=j
ite(nt)=j
maxite=max0(maxite,j)
if(maxite.gt.ntot)call bomb(4,'ntot ')
return
end
subroutine numl(j1,j2,j3,j4,j5)
c search j1-j2 find name in j3-j4, j5 is index of label
$include: 'pfcom.for'
character*1 c
if(j2.lt.j1)call bomb(3,'numl ')
call nxtch(j1,j2,j3)
do 10 j=j3,j2
c=a(j)
do 5 k=1,26
if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
5 continue
do 6 k=1,10
if(c.eq.hn(k))goto 10
6 continue
if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
10 continue
j4=j2
goto 14
12 j4=j-1
14 l=j4-j3+1
do 18 i=1,labels
if(ipl(i).ne.l)goto 18
do 16 j=1,l
if(pl(j,i).ne.a(j+j3-1))goto 18
16 continue
j5=i
return
18 continue
labels=labels+1
if(labels.gt.nlab)call bomb(4,'nlab ')
j5=labels
ipl(j5)=l
call copy(a(j3),pl(1,j5),l)
return
end
subroutine numv(j1,j2,j3,j4,j5)
c search j1-j2 find name in j3-j4, j5 is index of variable
$include: 'pfcom.for'
character*1 c
if(j2.lt.j1)call bomb(3,'numv ')
call nxtch(j1,j2,j3)
do 10 j=j3,j2
c=a(j)
do 5 k=1,26
if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
5 continue
do 6 k=1,10
if(c.eq.hn(k))goto 10
6 continue
if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
10 continue
j4=j2
goto 14
12 j4=j-1
14 l=j4-j3+1
do 28 ii=ns,1,-1
i=iosi(ii)
if(ipv(i).ne.l)goto 28
do 26 j=1,l
if(pv(j,i).ne.a(j+j3-1))goto 28
26 continue
j5=i
return
28 continue
do 18 i=1,nv
if(ittb(i).ne.0)goto 18
if(ipv(i).ne.l)goto 18
do 16 j=1,l
if(pv(j,i).ne.a(j+j3-1))goto 18
16 continue
j5=i
return
18 continue
nv=nv+1
j5=nv
if(nv.gt.nvar)call bomb(4,'nvar ')
ipv(nv)=l
call copy(a(j3),pv(1,nv),l)
return
end
subroutine nxt(j1,j2,h,j3)
character*1 h
call srch(j1,j2,h,j3)
if(j3.ne.-1)return
write(8,1)h,j1,j2
1 format(' nxt did not find ',a1,' in ',2i4)
call bomb(3,'nxt ')
return
end
subroutine nxtch(j1,j2,j3)
call srchch(j1,j2,j3)
if(j3.eq.-1)call bomb(3,'nxtch ')
return
end
subroutine nxtnb(j1,j2,j3)
call srchnb(j1,j2,j3)
if(j3.eq.-1)call bomb(3,'nxtnb ')
end
subroutine order(is)
$include: 'pfcom.for'
ina=ittb(is)
inb=itte(is)
if(inb.lt.ina .or. ina.eq.0)call bomb(3,'order1 ')
it=ina
2 ia=itb(it)
if(dabs(tm(it)).le.1.d-11)goto 41
ib=ite(it)
3 if(ib.eq.ia)goto 25
i=ia
4 if(ivp(i).eq.0)goto 20
i=i+1
if(i.le.ib)goto 4
if(ia.eq.ib)goto 25
23 i=ia+1
5 if(ivn(i)-ivn(i-1))10,7,6
6 i=i+1
if(i.le.ib)goto 5
goto 40
7 ivp(i-1)=ivp(i-1)+ivp(i)
if(i.eq.ib)goto 9
ib=ib-1
do 8 j=i,ib
ivn(j)=ivn(j+1)
8 ivp(j)=ivp(j+1)
if(ivp(i-1).ne.0)goto 5
goto 3
9 ib=ib-1
if(ivp(ib).ne.0)goto 40
goto 3
10 if(i.eq.ia+1)goto 19
ic=i-2
11 if(ivn(i)-ivn(ic))12,13,16
12 if(ic.eq.ia)goto 18
ic=ic-1
goto 11
13 ivp(ic)=ivp(ic)+ivp(i)
if(i.eq.ib)goto 15
ib=ib-1
do 14 j=i,ib
ivn(j)=ivn(j+1)
14 ivp(j)=ivp(j+1)
if(ivp(ic).eq.0)goto 3
goto 5
15 ib=ib-1
if(ivp(ic).ne.0)goto 40
goto 3
16 ip=ivp(i)
in=ivn(i)
do 17 j=i,ic+2,-1
ivn(j)=ivn(j-1)
17 ivp(j)=ivp(j-1)
ivn(ic+1)=in
ivp(ic+1)=ip
goto 6
18 ic=ia-1
goto 16
19 ip=ivp(i)
in=ivn(i)
ivp(i)=ivp(ia)
ivn(i)=ivn(ia)
ivp(ia)=ip
ivn(ia)=in
goto 6
20 if(i.eq.ib)goto 22
ib=ib-1
do 21 j=i,ib
ivn(j)=ivn(j+1)
21 ivp(j)=ivp(j+1)
goto 4
22 ib=ib-1
if(ib.gt.ia)goto 23
ib=ia
25 if(ivp(ia).eq.0)ivn(ia)=2
40 ite(it)=ib
if(ivn(ia).ne.2 .or. ivp(ia).le.maxe)goto 52
41 call bkup(it,inb)
goto 53
52 it=it+1
53 if(it.le.inb)goto 2
c sort terms ina-inb
if(ina.lt.inb)goto 70
if(ina.eq.inb)goto 82
51 inb=ina
tm(ina)=0.
k=itb(ina)
ite(ina)=k
ivn(k)=2
ivp(k)=0
goto 82
70 it=ina
68 it=it+1
69 if(it.gt.inb)goto 80
if(dabs(tm(it)).le.1.d-11)goto 75
if(icomp(ina,it))71,73,74
71 call tzap(it,ina)
goto 68
73 tm(ina)=tm(ina)+tm(it)
75 call bkup(it,inb)
goto 69
74 if(icomp(it-1,it))62,61,68
61 tm(it-1)=tm(it-1)+tm(it)
call bkup(it,inb)
goto 69
62 la=ina
lb=it-1
63 if(lb-la.eq.1)goto 64
lc=(la+lb)/2
if(icomp(lc,it))65,66,67
64 call tzap(it,lb)
goto 68
65 lb=lc
goto 63
66 tm(lc)=tm(lc)+tm(it)
call bkup(it,inb)
goto 69
67 la=lc
goto 63
80 if(ina.gt.inb)goto 51
82 call defsum(is,ina,inb)
return
end
subroutine ordset
$include: 'pfcom.for'
call nxt(2,nchpl,' ',j1)
call srchnb(j1,nchpl,j2)
if(j2.eq.-1)goto 20
if(a(j2).eq.';')goto 20
ntsav=nt
itesav=itetop
call term(j1)
nt=ntsav
itetop=itesav
return
20 if=1
do 10 il=1,300
call read(ib)
if(np.ge.2)call prntap
if(ib.eq.1)return
is=1
2 call srchch(is,nchpl,j1)
if(j1.eq.-1)goto 10
call numv(j1,nchpl,j2,j3,j4)
if=if+1
if(j4.lt.if)call bomb(3,'ordset1 ')
if(j4.eq.if)goto 8
do 5 i=nv,if,-1
ipv(i)=ipv(i-1)
do 4 ic=1,nch
4 pv(ic,i)=pv(ic,i-1)
5 continue
call copy(a(j2),pv(1,if),j3-j2+1)
ipv(if)=j3-j2+1
8 is=j3+1
if(is.le.nchpl)goto 2
10 continue
return
end
subroutine outqu(j5)
$include: 'pfcom.for'
1 do 6 it=ittb(j5),itte(j5)
do 5 iv=itb(it),ite(it)
if(ittb(ivn(iv)).ne.0)goto 9
5 continue
6 continue
return
9 jv=ivn(iv)
call subs(jv,jv,j5)
goto 1
end
subroutine pack
$include: 'pfcom.for'
if(minns.le.ns)goto 9
is=iosi(ns)
nt=itte(is)
itetop=ite(nt)
minns=ns+1
return
9 if(minns.le.0)call bomb(3,'pack4 ')
js=minns
kt=0
kv=0
if(js.eq.1)goto 13
kt=itte(iosi(js-1))
kv=ite(kt)
13 do 30 ms=js,ns
is=iosi(ms)
if(is.le.0)call bomb(3,'pack6 ')
ia=ittb(is)
ib=itte(is)
lv=0
lt=0
do 14 it=ia,ib
if(dabs(tm(it)).lt.1.d-11)goto 14
lt=lt+1
if(lt.gt.ntermt)call bomb(4,'ntermt ')
tmt(lt)=tm(it)
itbt(lt)=lv+1
do 12 i=itb(it),ite(it)
if(ivp(i).ne.0)goto 31
do 32 j=itb(it),ite(it)
if(ivp(j).ne.0)goto 12
32 continue
lv=lv+1
ifnt(lv)=2
ifpt(lv)=0
goto 33
31 lv=lv+1
ifnt(lv)=ivn(i)
ifpt(lv)=ivp(i)
12 continue
33 itet(lt)=lv
if(lv.gt.ntott)call bomb(4,'ntott ')
14 continue
if(lv.ne.0)goto 15
kt=kt+1
ittb(is)=kt
itte(is)=kt
kv=kv+1
tm(kt)=0.d0
itb(kt)=kv
ite(kt)=kv
ivn(kv)=2
ivp(kv)=0
goto 30
15 ittb(is)=kt+1
do 18 it=1,lt
kt=kt+1
tm(kt)=tmt(it)
itb(kt)=kv+1
do 16 i=itbt(it),itet(it)
kv=kv+1
ivn(kv)=ifnt(i)
16 ivp(kv)=ifpt(i)
18 ite(kt)=kv
itte(is)=kt
30 continue
nt=kt
itetop=kv
minns=ns+1
return
end
subroutine print(ii)
$include: 'pfcom.for'
common/pr1/pa(232),aa(22)
character*1 aa,pa,bl,sm,cf
data bl,sm,cf/' ','-','^'/
if(np.eq.0)return
ia=ittb(ii)
if(ia.eq.0)write(8,2)
if(ia.eq.0)return
2 format('there is nothing to print')
ib=itte(ii)
nst=21
do 40 it=ia,ib
n=nst
do 4 i=1,232
4 pa(i)=' '
ka=itb(it)
kb=ite(it)
do 30 k=ka,kb
n=n+1
pa(n)=bl
if=ivn(k)
jp=ivp(k)
ifl=ipv(if)
do 20 j=1,ifl
n=n+1
20 pa(n)=pv(j,if)
if(jp.eq.1)goto 30
n=n+1
if(jp.ge.0)pa(n)=cf
if(jp.lt.0)pa(n)=sm
iap=iabs(jp)
n=n+1
j=mod(iap,10)
pa(n)=hn(j+1)
if(iap.le.9)goto 30
j=iap/10
if(j.ge.10)call bomb(3,'print2 ')
n=n+1
pa(n)=pa(n-1)
pa(n-1)=hn(j+1)
30 continue
if(n.gt.232)call bomb(3,'print3 ')
call fmtc(it,if,il,' ')
if(il-if.eq.2 .and. aa(if+1).eq.'1' .and. aa(il).eq.'.')il=if
do 32 i=if,il
32 pa(nst-il+i)=aa(i)
j=it+1-ia
33 if(j.ge.1000)j=j-1000
if(j.ge.1000)goto 33
write(8,34)j,(pa(k),k=1,n)
34 format(i3,2x,127a1,/,27x,105a1)
40 continue
return
end
subroutine prntap
$include: 'pfcom.for'
do 1 n=nchpl,1,-1
if(ap(n).ne.' ')goto 2
1 continue
n=1
2 write(8,3)(ap(i),i=1,n)
3 format(131a1)
return
end
subroutine read(ib)
$include: 'pfcom.for'
c ib=0 on return if line has nonblanks
c ib=1 on return if line is all blanks
line=line+1
ib=inpba(line)
ie=inpea(line)
il=ie-ib+1
do 1 i=1,nchpl
1 a(i)=' '
call copy(inpl(ib),a,il)
call copy(a,ap,nchpl)
i=1
if(a(1).eq.'*')goto 8
do 6 i=1,nchpl
if(a(i).eq.';')goto 8
6 continue
goto 14
8 do 10 j=i,nchpl
10 a(j)=' '
14 ib=0
do 16 i=1,nchpl
if(a(i).ne.' ')goto 18
16 continue
17 ib=1
return
18 if(i.eq.1)return
do 20 j=i,nchpl
a(j+1-i)=a(j)
20 a(j)=' '
return
end
subroutine readin
$include: 'pfcom.for'
do 1 i=1,nlab
1 linlab(i)=-1
line=0
lines=0
labels=0
iinpl=1
5 read(7,6)a
6 format(127a1)
do 7 k=nchpl,1,-1
if(a(k).ne.' ')goto 10
7 continue
k=1
10 lines=lines+1
lintot=lintot+1
if(lines.gt.ninl)call bomb(4,'ninl ')
if(k.eq.1)goto 15
call srchnb(1,k,m)
if(a(m).eq.'p' .and. a(m+1).eq.'r' .and. a(m+2).eq.'o')goto 11
if(a(m).ne.'l' .or. a(m+1).ne.'a' .or. a(m+2).ne.'b')goto 15
11 call nxt(m+2,nchpl,' ',j1)
call nxtch(j1,nchpl,j3)
call nxt(j3,nchpl,' ',j4)
call numl(j1,nchpl,j3,j4,j5)
if(linlab(j5).ne.-1)write(8,12)(a(i),i=j3,j4)
12 format('duplicate label: ',15a1)
if(linlab(j5).ne.-1)call bomb(3,'readin1 ')
linlab(labels)=lines
15 inpba(lines)=iinpl
inpea(lines)=iinpl+k-1
if(inpea(lines).gt.ninch)call bomb(4,'ninch ')
call copy(a,inpl(iinpl),k)
iinpl=iinpl+k
if(k.eq.1)goto 5
if(a(m).eq.'f' .and. a(m+1).eq.'l' .and. a(m+2).eq.'u')goto 16
if(a(m).ne.'e' .or. a(m+1).ne.'n' .or. a(m+2).ne.'d')goto 5
if(a(m+3).ne.' ' .or. a(m+4).ne.'i' .or. a(m+5).ne.'n')goto 5
16 write(8,20)lines,labels,iinpl
20 format('number of lines in input file=',i5,' labels='
+ ,i2,' characters=',i6)
return
end
subroutine recip(j13,j4)
$include: 'pfcom.for'
ia=ittb(j13)
iz=itte(j13)
k=itb(ia)
if(ivn(k).eq.2 .and. ivp(k).ne.0 .and. ia.ne.iz)then
write(8,2)
call bomb(3,'recip1 ')
2 format('this simple program cannot process the denominator.',
+ ' denominators should have 1 term without variable 1.')
endif
ja=nt+1
call cpynt(ia,ia,1.d0)
tm(ja)=1.d0/tm(ja)
do 4 i=itb(ja),ite(ja)
4 ivp(i)=-ivp(i)
if(iz.eq.ia)call defsum(j4,ja,nt)
if(iz.eq.ia)return
if(maxe.gt.2)write(8,6)
if(maxe.gt.2)call bomb(3,'recip2 ')
6 format('this simple program cannot divide if maxe is',
+ ' greater than 2')
ib=ia+1
k=itb(ib)
if(ivn(k).ne.2 .or. ivp(k).lt.1)write(8,2)
if(ivn(k).ne.2 .or. ivp(k).lt.1)call bomb(3,'recip3 ')
if(ivp(k).eq.2)goto 14
ic=ib
if(ib.eq.iz)goto 10
8 k=itb(ic+1)
if(ic.eq.iz .or. ivp(k).eq.2)goto 10
ic=ic+1
goto 8
10 jb=nt+1
call cpynt(ib,ic,1.d0)
jc=nt
if(ic.eq.iz)goto 12
11 jd=nt+1
call cpynt(ic+1,iz,1.d0)
je=nt
goto 16
12 jd=nt+1
je=jd
itb(jd)=itetop+1
call incnt
call nite(itb(jd))
k=itetop
ivn(k)=2
ivp(k)=2
tm(jd)=0.d0
goto 16
14 jb=nt+1
jc=jb
itb(jb)=itetop+1
call incnt
call nite(itb(jb))
k=itetop
ivn(k)=2
ivp(k)=1
tm(jb)=0.d0
ic=ib-1
goto 11
16 jf=nt+1
call mult(ja,ja,jb,jc)
jg=nt
if(maxe.eq.1)goto 18
jh=nt+1
call mult(jf,jg,jf,jg)
ji=nt
jj=nt+1
call mult(ja,ja,jd,je)
jk=nt
jl=nt+1
call cpynt(jf,jg,-1.d0)
call cpynt(jh,ji,1.d0)
call cpynt(jj,jk,-1.d0)
jm=nt
call cpynt(ja,ja,1.d0)
call mult(ja,ja,jl,jm)
nt1=jm+1
goto 24
18 do 19 i=jf,jg
19 tm(i)=-tm(i)
call cpynt(ja,ja,1.d0)
call mult(ja,ja,jf,jg)
nt1=jg+1
24 call defsum(j4,nt1,nt)
call order(j4)
return
end
subroutine swap(j4)
$include: 'pfcom.for'
if=ittb(j4)
ig=ios(j4)
if(if.ne.0)minns=min0(minns,ios(j4))
ios(j4)=ios(1)
iosi(ios(j4))=j4
ittb(j4)=ittb(1)
itte(j4)=itte(1)
if(if.eq.0)goto 5
ios(1)=ig
iosi(ig)=1
ittb(1)=if
itte(1)=if
ite(if)=itb(if)
call pack
return
5 call niltrm
ittb(1)=0
call defsum(1,nt,nt)
return
end
subroutine srch(j1,j2,h,j3)
$include: 'pfcom.for'
character*1 h
do 1 j3=j1,j2
if(a(j3).eq.h)return
1 continue
j3=-1
return
end
subroutine srchch(j1,j2,j3)
$include: 'pfcom.for'
do 2 j3=j1,j2
do 1 i=1,26
if(a(j3).eq.ha(i))return
if(a(j3).eq.hc(i))return
1 continue
2 continue
j3=-1
end
subroutine srchnb(j1,j2,j3)
$include: 'pfcom.for'
do 1 j3=j1,j2
if(a(j3).ne.' ')return
1 continue
j3=-1
return
end
subroutine sub
$include: 'pfcom.for'
character*3 aj4,aj8
call nxt(1,nchpl,' ',j1)
call nxtnb(j1,nchpl,j2)
call nxtch(j1,nchpl,j3)
if(j2.ne.j3)goto 10
call numv(j1,nchpl,j2,j3,j4)
aj4='sum'
if(ittb(j4).eq.0)aj4='var'
1 call nxt(j3+1,nchpl,'f',j5)
if(a(j5+1).ne.'o')call bomb(3,'readin2 ')
if(a(j5+2).ne.'r')call bomb(3,'readin4 ')
if(a(j5+3).ne.' ')call bomb(3,'readin6 ')
call numv(j5+3,nchpl,j6,j7,j8)
aj8='var'
if(ittb(j8).gt.0)aj8='sum'
call nxt(j7+1,nchpl,'i',j9)
if(a(j9+1).ne.'n')call bomb(3,'readin8 ')
if(a(j9+2).ne.' ')call bomb(3,'readin10')
call numv(j9+2,nchpl,j10,j11,j12)
minns=min0(minns,ios(j12))
if(aj4.eq.'num')goto 35
if(aj4.eq.'var')goto 30
if(aj8.eq.'sum')goto 50
call subs(j4,j8,j12)
return
10 call coef(j1,j4,c)
aj4='num'
j3=j4-1
goto 1
30 if(aj8.eq.'sum')goto 40
call subv(j4,j8,j12)
return
35 if(aj8.eq.'sum')call bomb(3,'readin12')
call subn(c,j8,j12)
return
40 call subr(j8,j4,j12)
return
50 a(71)='A'
a(72)='?'
a(73)='0'
call numv(71,73,j15,j16,j17)
call subr(j8,j17,j12)
call subs(j4,j17,j12)
return
end
subroutine subn(c,j8,j12)
$include: 'pfcom.for'
if(ittb(j12).eq.0)call bomb(3,'subn1 ')
if(ittb(j8).ne.0)call bomb(3,'subn2 ')
do 8 it=ittb(j12),itte(j12)
do 5 i=itb(it),ite(it)
if(ivn(i).ne.j8)goto 5
tm(it)=tm(it)*c**ivp(i)
ivp(i)=0
if(ite(it).eq.itb(it))goto 5
ite(it)=ite(it)-1
if(i.gt.ite(it))goto 5
do 1 j=i,ite(it)
ivn(j)=ivn(j+1)
1 ivp(j)=ivp(j+1)
5 continue
8 continue
call order(j12)
return
end
subroutine subr(ia,ib,ic)
c replace pattern=first term in sum ia by variable ib in sum ic
$include: 'pfcom.for'
if(ittb(ia).eq.0)call bomb(3,'subr1 ')
if(ittb(ib).ne.0)call bomb(3,'subr2 ')
nt1=nt+1
j1=ittb(ia)
j2=itb(j1)
j3=ite(j1)
j4=ittb(ic)
j5=itte(ic)
do 30 it=j4,j5
j6=itb(it)
j7=ite(it)
minrat=9999
do 5 j8=j2,j3
jfp=ivn(j8)
1 jfm=ivn(j6)
if(jfm.gt.jfp)goto 22
irat=ivp(j6)/ivp(j8)
if(jfm.eq.jfp .and. irat.ge.1)goto 4
jt=it+1-j4
if(j2.eq.j3 .and. np.ge.2 .and. jfm.eq.jfp .and. irat.le.-1)
+ write(8,2)jt,(pv(i,ic),i=1,nch),(pv(i,jfp),i=1,nch),
+ (pv(i,ia),i=1,nch)
2 format('term',i4,' of ',15a1,' has exponent of ',15a1,
+ ' with different sign from ',15a1)
j6=j6+1
if(j6.le.j7)goto 1
goto 22
4 minrat=min0(minrat,ivp(j6)/ivp(j8))
5 continue
call incnt
tm(nt)=tm(it)/tm(j1)
itb(nt)=itetop+1
j6=itb(it)
k=itb(nt)
do 20 j=j6,j7
do 12 j8=j2,j3
if(ivn(j8).eq.ivn(j))goto 16
12 continue
ivn(k)=ivn(j)
ivp(k)=ivp(j)
goto 20
16 ivn(k)=ivn(j)
ivp(k)=ivp(j)-minrat*ivp(j8)
20 k=k+1
ivn(k)=ib
ivp(k)=minrat
call nite(k)
goto 30
22 call cpynt(it,it,1.d0)
30 continue
call defsum(ic,nt1,nt)
call order(ic)
return
end
subroutine subv(j4,j8,j12)
$include: 'pfcom.for'
if(ittb(j12).eq.0)call bomb(3,'subv1 ')
if(ittb(j4).ne.0 .or. ittb(j8).ne.0)call bomb(3,'subv2 ')
do 8 it=ittb(j12),itte(j12)
do 5 i=itb(it),ite(it)
if(ivn(i).eq.j8)ivn(i)=j4
5 continue
8 continue
call order(j12)
return
end
subroutine subs(in1,in2,in3)
c substitute sum=in1 for variable=in2 in sum=in3
c result goes in nt+1,...
$include: 'pfcom.for'
common/nwt/iitb,iite,jitb,jite,ntb,if
c dimension of ittbs ittes is max jp
dimension ittbs(30),ittes(30)
jp=1
ia=in2
ib=ittb(in3)
ic=itte(in3)
id=ittb(in1)
ie=itte(in1)
if(ib.eq.0 .or. id.eq.0)call bomb(3,'subs1 ')
ntb=nt+1
lev1=0
do 30 it=ib,ic
iitb=itb(it)
iite=ite(it)
if(ivn(iitb).eq.2)lev1=ivp(iitb)
if(lev1.gt.maxe)goto 31
do 5 if=iitb,iite
if(ivn(if).ne.ia)goto 5
ifp=ivp(if)
lev2=0
if(ifp.eq.1)goto 10
if(id.eq.ie)goto 70
if(ifp.gt.0)goto 17
itpr=it-ib+1
if(np.ge.2)write(8,1)(pv(i,ia),i=1,nch),
+ (pv(i,in3),i=1,nch),itpr
1 format(2x,15a1,' occurs to a negative power in ',15a1,
+ ' term ',i4)
5 continue
c here copy term it to nt+1
if=0
jitb=0
call sub2(tm(it))
goto 30
c here copy term it with substitution
10 do 68 jt=id,ie
jitb=itb(jt)
jite=ite(jt)
if(ivn(jitb).eq.2)lev2=ivp(jitb)
if(lev1+lev2.gt.maxe)goto 30
call sub2(tm(it)*tm(jt))
68 continue
goto 30
70 if(ivn(itb(id)).eq.2)lev2=ivp(itb(id))*ifp
if(lev1+lev2.gt.maxe)goto 30
k=0
do 71 i=itb(id),ite(id)
k=k+1
ifnt(k)=ivn(i)
71 ifpt(k)=ivp(i)*ifp
jitb=1
jite=k
call sub3(tm(it)*tm(id)**ifp)
goto 30
c here get in1**ifp
17 jd=id
je=ie
ifpsav=ifp
ntbsav=ntb
iitbsv=iitb
iitesv=iite
ifsav=if
ntbmu=nt
itesav=itetop
c test whether this power is in memory
if(jp.eq.1)goto 80
if(ifp.lt.jp)goto 27
if(ittbs(jp).eq.0)goto 30
if(ifp.eq.jp)goto 27
jd=nt+1
k=itetop
do 86 i=ittbs(jp),ittes(jp)
call incnt
tm(nt)=tmt(i)
itb(nt)=k+1
do 84 j=itbt(i),itet(i)
k=k+1
ivn(k)=ifnt(j)
84 ivp(k)=ifpt(j)
call nite(k)
86 continue
je=nt
goto 18
80 ks=0
nts=0
18 jpb=nt+1
ntbtmu=nt
call mult(id,ie,jd,je)
if(nt.eq.ntbtmu)then
jp=jp+1
ittbs(jp)=0
nt=ntbmu
ntb=ntbsav
itetop=itesav
goto 30
endif
jpe=nt
jp=jp+1
if(jp.gt.30)call bomb(4,'max jp ')
jd=jpb
je=jpe
ittbs(jp)=nts+1
do 21 i=jpb,jpe
nts=nts+1
tmt(nts)=tm(i)
itbt(nts)=ks+1
do 20 j=itb(i),ite(i)
ks=ks+1
ifnt(ks)=ivn(j)
20 ifpt(ks)=ivp(j)
itet(nts)=ks
if(ks.gt.ntott)call bomb(4,'ntott ')
if(nts.gt.ntermt)call bomb(4,'ntermt ')
21 continue
ittes(jp)=nts
if(jp.ne.ifp)goto 18
nt=ntbmu
itetop=itesav
if=ifsav
ifp=ifpsav
iitb=iitbsv
iite=iitesv
ntb=ntbsav
27 do 54 jt=ittbs(ifp),ittes(ifp)
jitb=itbt(jt)
jite=itet(jt)
if(ifnt(jitb).eq.2)lev2=ifpt(jitb)
if(lev1+lev2.gt.maxe)goto 30
call sub3(tm(it)*tmt(jt))
54 continue
30 continue
31 call defsum(in3,ntb,nt)
return
end
subroutine sub2(coef)
$include: 'pfcom.for'
common/nwt/iitb,iite,jitb,jite,ntb,if
if(dabs(coef).lt.1.d-11)return
call incnt
k=itetop
itb(nt)=k+1
tm(nt)=coef
ii=iitb
if(jitb.eq.0)goto 15
jj=jitb
11 if(ii.gt.iite)goto 13
if(ii.eq.if)goto 12
if(jj.gt.jite)goto 15
if(ivn(ii).lt.ivn(jj))goto 16
if(ivn(ii).eq.ivn(jj))goto 57
if(ivp(jj).eq.0)goto 18
k=k+1
ivn(k)=ivn(jj)
ivp(k)=ivp(jj)
18 jj=jj+1
goto 11
12 ii=ii+1
goto 11
13 if(jj.gt.jite)goto 60
if(ivp(jj).eq.0)goto 17
k=k+1
ivn(k)=ivn(jj)
ivp(k)=ivp(jj)
17 jj=jj+1
goto 13
14 if(ii.gt.iite)goto 60
15 if(ivp(ii).eq.0)goto 9
k=k+1
ivn(k)=ivn(ii)
ivp(k)=ivp(ii)
9 ii=ii+1
if(ii.eq.if)goto 9
goto 14
16 if(ivp(ii).eq.0)goto 12
k=k+1
ivn(k)=ivn(ii)
ivp(k)=ivp(ii)
ii=ii+1
goto 11
57 if(ivp(ii)+ivp(jj).eq.0)goto 58
k=k+1
ivn(k)=ivn(ii)
ivp(k)=ivp(ii)+ivp(jj)
58 ii=ii+1
jj=jj+1
goto 11
60 if(k.gt.itetop)goto 56
k=itetop+1
ivn(k)=2
ivp(k)=0
56 itesav=itetop
call nite(k)
if(nt.eq.ntb)goto 68
if(icomp(ntb,nt))31,33,34
31 call tzap(nt,ntb)
goto 68
33 tm(ntb)=tm(ntb)+tm(nt)
nt=nt-1
itetop=itesav
goto 68
34 if(icomp(nt-1,nt))62,61,68
61 tm(nt-1)=tm(nt-1)+tm(nt)
nt=nt-1
itetop=itesav
goto 68
62 la=ntb
lb=nt-1
63 if(lb-la.eq.1)goto 64
lc=(la+lb)/2
if(icomp(lc,nt))65,66,67
64 call tzap(nt,lb)
goto 68
65 lb=lc
goto 63
66 tm(lc)=tm(lc)+tm(nt)
nt=nt-1
itetop=itesav
goto 68
67 la=lc
goto 63
68 return
end
subroutine sub3(coef)
$include: 'pfcom.for'
common/nwt/iitb,iite,jitb,jite,ntb,if
if(dabs(coef).lt.1.d-11)return
call incnt
k=itetop
itb(nt)=k+1
tm(nt)=coef
ii=iitb
if(jitb.eq.0)goto 15
jj=jitb
11 if(ii.gt.iite)goto 13
if(ii.eq.if)goto 12
if(jj.gt.jite)goto 15
if(ivn(ii).lt.ifnt(jj))goto 16
if(ivn(ii).eq.ifnt(jj))goto 57
if(ifpt(jj).eq.0)goto 18
k=k+1
ivn(k)=ifnt(jj)
ivp(k)=ifpt(jj)
18 jj=jj+1
goto 11
12 ii=ii+1
goto 11
13 if(jj.gt.jite)goto 60
if(ifpt(jj).eq.0)goto 17
k=k+1
ivn(k)=ifnt(jj)
ivp(k)=ifpt(jj)
17 jj=jj+1
goto 13
14 if(ii.gt.iite)goto 60
15 if(ivp(ii).eq.0)goto 9
k=k+1
ivn(k)=ivn(ii)
ivp(k)=ivp(ii)
9 ii=ii+1
if(ii.eq.if)goto 9
goto 14
16 if(ivp(ii).eq.0)goto 12
k=k+1
ivn(k)=ivn(ii)
ivp(k)=ivp(ii)
ii=ii+1
goto 11
57 if(ivp(ii)+ifpt(jj).eq.0)goto 58
k=k+1
ivn(k)=ivn(ii)
ivp(k)=ivp(ii)+ifpt(jj)
58 ii=ii+1
jj=jj+1
goto 11
60 if(k.gt.itetop)goto 56
k=itetop+1
ivn(k)=2
ivp(k)=0
56 itesav=itetop
call nite(k)
if(nt.eq.ntb)goto 68
if(icomp(ntb,nt))31,33,34
31 call tzap(nt,ntb)
goto 68
33 tm(ntb)=tm(ntb)+tm(nt)
nt=nt-1
itetop=itesav
goto 68
34 if(icomp(nt-1,nt))62,61,68
61 tm(nt-1)=tm(nt-1)+tm(nt)
nt=nt-1
itetop=itesav
goto 68
62 la=ntb
lb=nt-1
63 if(lb-la.eq.1)goto 64
lc=(la+lb)/2
if(icomp(lc,nt))65,66,67
64 call tzap(nt,lb)
goto 68
65 lb=lc
goto 63
66 tm(lc)=tm(lc)+tm(nt)
nt=nt-1
itetop=itesav
goto 68
67 la=lc
goto 63
68 return
end
subroutine switch(j,k,is)
$include: 'pfcom.for'
do 2 it=ittb(is),itte(is)
do 1 i=itb(it),ite(it)
if(ivn(i).eq.j)ivn(i)=k
1 continue
2 continue
return
end
subroutine taylor
$include: 'pfcom.for'
character*4 star
dimension ltay(mitay)
data iorder,ifirst/1,0/
if(idflg.ne.0)call dertab
call nxt(2,nchpl,' ',j1)
call srch(j1,nchpl,'=',j2)
j4=j1
if(j2.eq.-1)goto 5
call getexp(j2+1,iorder,j3,0)
if(iorder.lt.1)call bomb(3,' taylor1')
return
5 call nxt(2,nchpl,' ',j13)
star='no'
if(a(j13-1).eq.'*')star='yes'
10 call nxt(j4,nchpl,'i',j3)
if(j3.eq.-1)call bomb(3,' taylor2')
if(a(j3-1).eq.' ' .and. a(j3+1).eq.'n' .and. a(j3+2).eq.' ')
+ goto 12
j4=j3+1
if(j4.lt.nchpl)goto 10
call bomb(3,' taylor3')
12 do 14 j5=j3-1,j1,-1
if(a(j5).ne.' ')goto 16
14 continue
call bomb(3,' taylor4')
16 do 18 j6=j5,j1,-1
if(a(j6).eq.' ')goto 20
18 continue
call bomb(3,' taylor5')
20 call numv(j6,j5,j7,j8,j9)
if(ittb(j9).eq.0)call bomb(3,' taylor6')
itesav=itetop
call term(j3+2)
itay=itetop-itesav
if(itay.gt.mitay)call bomb(4,'mitay ')
do 21 i=1,itay
21 ltay(i)=ivn(itesav+i)
nt=nt-1
itetop=itesav
call niltrm
call defsum(1,nt,nt)
if(ifirst.ne.0)goto 22
a(1)='A'
a(2)='?'
a(3)='0'
call numv(1,3,j14,j15,j10)
a(3)='1'
call numv(1,3,j14,j15,j11)
a(3)='2'
call numv(1,3,j14,j15,j12)
22 ifirst=1
nv1=nv+1
do 34 it=ittb(j9),itte(j9)
call cpynt(it,it,1.d0)
call defsum(j10,nt,nt)
do 32 ivv=1,itay
iv=ltay(ivv)
if(star.eq.'yes')call switch(iv,nv1,j10)
call copys(j10,j11)
call subn(0.d0,iv,j10)
f=1.d0
do 30 io=1,iorder
call deriv(j11,iv,j12)
call outqu(j12)
if(tm(ittb(j12)).eq.0.d0)goto 31
call copys(j12,j11)
call subn(0.d0,iv,j12)
f=f*io
call incnt
k=itetop+1
tm(nt)=1.d0/f
itb(nt)=k
ite(nt)=k
ivn(k)=iv
ivp(k)=io
call nite(k)
nt1=nt+1
call mult(nt,nt,ittb(j12),itte(j12))
call defsum(j12,nt1,nt)
call pack
30 call add(j12,j10,1.d0,j10)
31 if(star.eq.'yes')call switch(nv1,iv,j10)
32 continue
34 call add(j10,1,1.d0,1)
call swap(j9)
call order(j9)
call delsum(j10)
return
end
subroutine term(ist)
c each term is on a separate line, free format with coefficient first
$include: 'pfcom.for'
call incnt
call coef(ist,j1,c)
tm(nt)=c
itb(nt)=itetop+1
if(j1.eq.-1)then
k=itb(nt)
ivn(k)=2
ivp(k)=0
goto 10
endif
k=itetop
5 ise=1
if(a(j1).eq.'*')j1=j1+1
if(a(j1).eq.'/')ise=-1
if(a(j1).eq.'/')j1=j1+1
call numv(j1,nchpl,j2,j3,j4)
ie=1
j5=-1
if(j3.eq.nchpl-1 .and. a(nchpl).ne.' ')call bomb(3,'term1 ')
if(j3.le.nchpl-2)call getexp(j3+1,ie,j5,1)
k=k+1
ivn(k)=j4
ivp(k)=ie*ise
j1=j5
if(j5.ne.-1)goto 5
10 call nite(k)
return
end
subroutine thderv
$include: 'pfcom.for'
idflg=1
call nxt(2,nchpl,'o',j1)
call nxt(j1,nchpl,' ',j2)
call numv(j2,nchpl,j3,j4,j5)
call nxt(j4+1,nchpl,'w',j6)
call nxt(j6,nchpl,' ',j7)
call numv(j7,nchpl,j8,j9,j10)
call srch(j9+1,nchpl,'=',j11)
if(j11.ne.-1)goto 1
call srch(j9+1,nchpl,'i',j12)
if(j12.eq.-1)call bomb(3,'thderv1 ')
call nxt(j12,nchpl,' ',j11)
1 call srchnb(j11+1,nchpl,j12)
if(j12.eq.-1)call bomb(3,'thderv2 ')
ig=0
if(a(j12).eq.'0')ig=1
if(ig.eq.1)goto 8
call numv(j11+1,nchpl,j12,j13,j14)
8 if=0
do 10 i=1,nd
if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)if=1
if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)in=i
10 continue
if(if.eq.0 .and. ig.eq.1)goto 50
if(if.eq.0)goto 14
ja=ider(3,in)
jb=ipv(ja)
if(np.ge.2)write(8,11)(pv(i,ja),i=1,jb)
11 format('REDEFINITION: previously the derivative was = ',15a1)
if(ig.eq.0)goto 18
nd=nd-1
if(in.gt.nd)goto 50
do 28 ih=in,nd
ider(1,ih)=ider(1,ih+1)
ider(2,ih)=ider(2,ih+1)
28 ider(3,ih)=ider(3,ih+1)
goto 50
14 nd=nd+1
in=nd
if(nd.gt.nder)call bomb(4,'nder ')
ider(1,in)=j5
ider(2,in)=j10
18 ider(3,in)=j14
50 return
end
subroutine tzap(ka,kb)
c zap term at ka to position kb
$include: 'pfcom.for'
if(ka.eq.kb)return
tmka=tm(ka)
itbka=itb(ka)
iteka=ite(ka)
if(ka.gt.kb)goto 5
do 1 i=ka,kb-1
tm(i)=tm(i+1)
itb(i)=itb(i+1)
1 ite(i)=ite(i+1)
goto 8
5 do 6 i=ka,kb+1,-1
tm(i)=tm(i-1)
itb(i)=itb(i-1)
6 ite(i)=ite(i-1)
8 tm(kb)=tmka
itb(kb)=itbka
ite(kb)=iteka
return
end
subroutine zsub(j3)
return
end